function [fval, c_star, K, Y, G, Z, H, K_bar, E_a, E_p, q_capital, q_money, q_profit, w, ...
    L_s, L_b, ALPHA_0, ALPHA_0_meet, ALPHA_0_tilde_meet, Pss_tilde_meet, K_meet, K_tilde_meet, a_meet, a_tilde_meet,...
    int_aqc, int_sppe, int_prob_aqc, int_prob_target, int_average_2nd_price, int_prod_std, int_HT_credit, int_debt, q_surplus_entry, int_surplus_buy, int_surplus_sell] =...
    eqm_ss_persistent_case_6_log_normal(x, par, fun, eqm_all)

%%%% check if endogenous entry is shut down
%%%% functions always have seller' productivity in the 1st argument
%%%% cutoff value indicates when the division line reaches the upper bound of buyer's productivity

int_method  = par.int_method;
Rel_Tol     = par.Rel_Tol;
Abs_Tol     = par.Abs_Tol;

n_s         = 2;            % number of states in the Markov matrix
B           = x(end);       % we always put B as the last component, before parameters to be calibrated       
Z           = zeros(n_s, 1);
K           = zeros(n_s, 1);
ALPHA_0     = zeros(n_s, 1);
fval        = zeros(n_s  * 3 + 1, 1);

for ii = 1:n_s
    Z(ii) = x(ii); 
    K(ii) = x(ii + n_s);
    ALPHA_0(ii) = x(ii + 2 * n_s);
end
w       = par.ETA / (B / par.A / (1 - par.ETA))^((1 - par.ETA) / par.ETA);


CHI_q       = par.CHI_q;
CHI_PI      = par.CHI_PI;
CHI_PI_k    = par.CHI_PI_k;
CHI_k       = par.CHI_k;
XI          = par.XI;
ALPHA       = par.ALPHA;
G           = par.G;

Z_meet          = Z(par.state(:,1)); % the first half is high state; the second half is low state
Z_tilde_meet    = Z(par.state_tilde(:,1)); % the repetition is high high
K_meet          = K(par.state(:,1)); % the first half is high state; the second half is low state
K_tilde_meet    = K(par.state_tilde(:,1));
a_meet          = par.a(par.state(:,2));
a_tilde_meet    = par.a(par.state_tilde(:,2));

ALPHA_0_meet        = ALPHA_0(par.state(:,1));
ALPHA_0_tilde_meet  = ALPHA_0(par.state_tilde(:,1));
Pss_tilde_meet      = par.Pss_tilde_meet.* ALPHA_0_tilde_meet / (par.Pss' * ALPHA_0);


%%%% log(a) is the constant used to adjust the log normal
%%%% distribution, the lower bound and upper bound stay the same
%%%% note: unlike the iid case, we do not solve for L directly; so we could
%%%% bound L in the following (but not in the integrals below)
L_s             = max(min(((Z_tilde_meet + CHI_k * (1 - par.DELTA) * K_tilde_meet) ./ K_meet - (1 - par.DELTA) * (1 - CHI_q))...
                    ./ (B * (1 - par.tau_k)), par.EPS_H_BOUND), par.EPS_L_BOUND);
L_b             = max(min(((Z_meet + CHI_k * (1 - par.DELTA) * K_meet) ./ K_tilde_meet - (1 - par.DELTA) * (1 - CHI_q))...
                    ./ (B * (1 - par.tau_k)), par.EPS_H_BOUND), par.EPS_L_BOUND);


%%%% Considering the vector of integrals
int_capital_gain= arrayfun(@(L, a_s, a_b, K_K_tilde) integral2(@(eps, eps_t)fun.capital_gain(eps, eps_t, a_s, a_b),...
    fun.l_bound(a_s), min(max(L / (1 - CHI_PI - K_K_tilde * CHI_PI_k), fun.l_bound(a_s)), fun.h_bound(a_s)),...
    @(eps)eps, @(eps)max(eps, fun.threshold(eps, L, CHI_PI, CHI_PI_k, K_K_tilde)),...
    'method', int_method, 'RelTol', Rel_Tol, 'AbsTol', Abs_Tol),...
    L_s, a_meet, a_tilde_meet, K_tilde_meet./K_meet); %note for K_K_tilde, it means buyers' capital / sellers' capital

int_money_gain  = arrayfun(@(L, a_s, a_b, K_K_tilde) integral2(@(eps_t, eps)fun.money_gain(eps_t, eps, a_s, a_b, B, CHI_q, CHI_PI),...
    min(max(L / par.THETA - (1 - par.THETA - (1 + K_K_tilde) * CHI_PI) / par.THETA * fun.h_bound(a_b), fun.l_bound(a_s)), fun.h_bound(a_s)),...
    min(max(L / (1 - (1 + K_K_tilde) * CHI_PI), fun.l_bound(a_s)), fun.h_bound(a_s)),...
    @(eps_t)fun.threshold(eps_t, L, CHI_PI, CHI_PI_k, K_K_tilde), @(eps_t)eps_t - eps_t + fun.h_bound(a_b),...
    'method', int_method, 'RelTol', Rel_Tol, 'AbsTol', Abs_Tol)...
    + integral2(@(eps_t, eps) fun.money_gain(eps_t, eps, a_s, a_b, B, CHI_q, CHI_PI),...
    min(max(L / (1 - (1 + K_K_tilde) * CHI_PI), fun.l_bound(a_s)), fun.h_bound(a_s)), fun.h_bound(a_s),...
    @(eps_t)eps_t, @(eps_t)eps_t - eps_t + fun.h_bound(a_b),...
    'method',int_method,'RelTol', Rel_Tol, 'AbsTol', Abs_Tol),...
    L_b, a_tilde_meet, a_meet, K_meet./K_tilde_meet); %note for K_K_tilde, it means buyers' capital / sellers' capital

%%%% for Holmostrom - Tirol benefit
int_profit_gain  = arrayfun(@(L, a_s, a_b, K_K_tilde) integral2(@(eps_t, eps)fun.profit_gain(eps_t, eps, a_s, a_b, B, CHI_q, CHI_PI),...
    min(max(L / par.THETA - (1 - par.THETA - (1 + K_K_tilde) * CHI_PI) / par.THETA * fun.h_bound(a_b), fun.l_bound(a_s)), fun.h_bound(a_s)),...
    min(max(L / (1 -  (1 + K_K_tilde) * CHI_PI), fun.l_bound(a_s)), fun.h_bound(a_s)),...
    @(eps_t)fun.threshold(eps_t, L, CHI_PI, CHI_PI_k, K_K_tilde), @(eps_t)eps_t - eps_t + fun.h_bound(a_b),...
    'method', int_method, 'RelTol', Rel_Tol, 'AbsTol', Abs_Tol)...
    + integral2(@(eps_t, eps) fun.profit_gain(eps_t, eps, a_s, a_b, B, CHI_q, CHI_PI),...
    min(max(L / (1 -  (1 + K_K_tilde) * CHI_PI), fun.l_bound(a_s)), fun.h_bound(a_s)), fun.h_bound(a_s),...
    @(eps_t)eps_t, @(eps_t)eps_t - eps_t + fun.h_bound(a_b),...
    'method',int_method,'RelTol', Rel_Tol, 'AbsTol', Abs_Tol),...
    L_b, a_tilde_meet, a_meet, K_meet./K_tilde_meet); %note for K_K_tilde, it means buyers' capital / sellers' capital

int_liq   = arrayfun(@(L, a_s, a_b, Z_K, K_K_tilde) integral2(@(eps_t, eps)fun.liq_gain(eps_t, eps, a_s, a_b, B, CHI_q, CHI_PI, CHI_PI_k, CHI_k, Z_K),...
    min(max(L / par.THETA - (1 - par.THETA - (1 + K_K_tilde) * CHI_PI) / par.THETA *  fun.h_bound(a_b), fun.l_bound(a_s)), fun.h_bound(a_s)),...
    min(max(L / (1 - (1 + K_K_tilde) * CHI_PI), fun.l_bound(a_s)), fun.h_bound(a_s)),...
    @(eps_t)fun.threshold(eps_t, L, CHI_PI, CHI_PI_k, K_K_tilde), @(eps_t)eps_t - eps_t + fun.h_bound(a_b),...
    'method', int_method, 'RelTol', Rel_Tol, 'AbsTol', Abs_Tol)...
    + integral2(@(eps_t,eps)fun.liq_gain(eps_t, eps, a_s, a_b, B, CHI_q, CHI_PI, CHI_PI_k, CHI_k, Z_K),...
    min(max(L / (1 - (1 + K_K_tilde) * CHI_PI), fun.l_bound(a_s)), fun.h_bound(a_s)), fun.h_bound(a_s),... 
    @(eps_t)eps_t, @(eps_t)eps_t - eps_t + fun.h_bound(a_b),...
    'method', int_method, 'RelTol', Rel_Tol, 'AbsTol', Abs_Tol),...
    L_b, a_tilde_meet, a_meet, Z_meet./K_meet, K_meet./K_tilde_meet); %note for K_K_tilde, it means buyers' capital / sellers' capital

int_sppe        = arrayfun(@(L, a_s, a_b, Z_K, K_K_tilde)integral2(@(eps_t,eps)fun.reall_P(eps_t, eps, a_s, a_b, B, CHI_q, CHI_PI, CHI_PI_k, CHI_k, Z_K),...
    max(L / par.THETA - (1 - par.THETA - (1 + K_K_tilde) * CHI_PI) / par.THETA *  fun.h_bound(a_b), fun.l_bound(a_s)),...
    min(max(L / (1 - (1 + K_K_tilde) * CHI_PI), fun.l_bound(a_s)), fun.h_bound(a_s)),...
    @(eps_t)fun.threshold(eps_t, L, CHI_PI, CHI_PI_k, K_K_tilde), @(eps_t)eps_t - eps_t + fun.h_bound(a_b),...
    'method',int_method,'RelTol', Rel_Tol, 'AbsTol', Abs_Tol)...
    + integral2(@(eps_t, eps)fun.reall_P(eps_t, eps, a_s, a_b, B, CHI_q, CHI_PI, CHI_PI_k, CHI_k, Z_K),...
    min(max(L / (1 - (1 + K_K_tilde) * CHI_PI), fun.l_bound(a_s)), fun.h_bound(a_s)), fun.h_bound(a_s),...
    @(eps_t)eps_t, @(eps_t)eps_t - eps_t + fun.h_bound(a_b),...
    'method',int_method,'RelTol', Rel_Tol, 'AbsTol', Abs_Tol),...
    L_b, a_tilde_meet, a_meet, Z_meet./K_meet, K_meet./K_tilde_meet);  %note for K_K_tilde, it means buyers' capital / sellers' capital    

int_aqc         = arrayfun(@(L, a_s, a_b, K_K_tilde) integral2(@(eps, eps_t)fun.reall_L(eps, eps_t, a_s, a_b, B, CHI_q, CHI_PI, CHI_PI_k),...
    fun.l_bound(a_s), min(max(L / (1 - (1 + K_K_tilde) * CHI_PI), fun.l_bound(a_s)), fun.h_bound(a_s)),...
    @(eps)eps, @(eps)max(eps, fun.threshold(eps, L, CHI_PI, CHI_PI_k, K_K_tilde)),...
    'method', int_method, 'RelTol', Rel_Tol, 'AbsTol', Abs_Tol),...
    L_s, a_meet, a_tilde_meet, K_tilde_meet./K_meet); %note for K_K_tilde, it means buyers' capital / sellers' capital

%%%% Considering firms from different states
q_capital  = [sum(par.P_tr_meet(1:end/n_s) .* par.P_tr_tilde_meet(1:end/n_s) .* Pss_tilde_meet(1:end/n_s) .* int_capital_gain(1:end/n_s));...
     sum(par.P_tr_meet(end/n_s+1:end) .* par.P_tr_tilde_meet(end/n_s+1:end) .* Pss_tilde_meet(end/n_s+1:end) .* int_capital_gain(end/n_s+1:end))];           

q_money    = [sum(par.P_tr_meet(1:end/n_s) .* par.P_tr_tilde_meet(1:end/n_s) .* Pss_tilde_meet(1:end/n_s) .* int_money_gain(1:end/n_s));...
    sum(par.P_tr_meet(end/n_s+1:end) .* par.P_tr_tilde_meet(end/n_s+1:end) .* Pss_tilde_meet(end/n_s+1:end) .* int_money_gain(end/n_s+1:end))];

q_profit   = [sum(par.P_tr_meet(1:end/n_s) .* par.P_tr_tilde_meet(1:end/n_s) .* Pss_tilde_meet(1:end/n_s) .* int_profit_gain(1:end/n_s));...
    sum(par.P_tr_meet(end/n_s+1:end) .* par.P_tr_tilde_meet(end/n_s+1:end) .* Pss_tilde_meet(end/n_s+1:end) .* int_profit_gain(end/n_s+1:end))];

%%%% the following matrix manipulation turns out to be slower
% q_capital     = yesterday * (par.P_tr_tilde_meet .* Pss_tilde_meet .* int_capital_gain); 
% q_money       = yesterday * (par.P_tr_tilde_meet .* Pss_tilde_meet .* int_money_gain);


E_a = sum(ALPHA_0_meet .* par.P_tr_meet .* par.P_tr_tilde_meet .* par.Pss_meet .* Pss_tilde_meet .* K_meet .* int_aqc)  * ALPHA;        

E_p = sum(ALPHA_0_meet .* par.P_tr_meet .* par.P_tr_tilde_meet .* par.Pss_meet .* Pss_tilde_meet .* K_meet .* int_sppe) * ALPHA;


term_1 = sum(par.Pss .* K .* (par.P * par.q_average_prod));
term_2 = ALPHA * ...
    sum(ALPHA_0_meet .* par.P_tr_meet .* par.P_tr_tilde_meet .* par.Pss_meet .* Pss_tilde_meet .* K_meet .* int_capital_gain);
term_3 = ALPHA * ...
    sum(ALPHA_0_meet .* par.P_tr_meet .* par.P_tr_tilde_meet .* par.Pss_meet .* Pss_tilde_meet .* K_meet .* int_liq);
K_bar   = term_1 +  term_2 + term_3;


c_star  = (XI / (1 - par.tau_h) / w)^(1 / (-par.SIGMA));

Y       = B * K_bar / (1 - par.ETA);

goods_mkt = c_star + G + par.DELTA * sum(par.Pss .* K) - Y; % the goods mkt clearing condition 
          
int_surplus_buy = par.THETA * XI / w * arrayfun(@(L, a_s, a_b, Z_K, K_tilde_K)...
                integral2(@(eps_t, eps)fun.surplus_buy(eps_t, eps, a_s, a_b, B, CHI_q, CHI_PI, CHI_PI_k, CHI_k, Z_K, K_tilde_K),...
                fun.l_bound(a_s), fun.h_bound(a_s), @(eps_t)eps_t, fun.h_bound(a_b),...
                'method', int_method, 'RelTol', Rel_Tol, 'AbsTol', Abs_Tol),...
                L_b, a_tilde_meet, a_meet, Z_meet./ K_meet, K_tilde_meet./K_meet) .* K_meet;

int_surplus_sell= (1 - par.THETA) * XI / w * arrayfun(@(L, a_s, a_b, Z_K, K_tilde_K)...
                integral2(@(eps, eps_t)fun.surplus_buy(eps, eps_t, a_s, a_b, B, CHI_q, CHI_PI, CHI_PI_k, CHI_k, Z_K, K_tilde_K),...
                fun.l_bound(a_s), fun.h_bound(a_s), @(eps)eps, fun.h_bound(a_b),...
                'method', int_method, 'RelTol', Rel_Tol, 'AbsTol', Abs_Tol),...
                L_s, a_meet, a_tilde_meet, Z_tilde_meet./K_tilde_meet, K_meet./K_tilde_meet) .* K_tilde_meet;

q_surplus_entry   = ...
      [sum(par.P_tr_meet(1:end/n_s) .* par.P_tr_tilde_meet(1:end/n_s) .* Pss_tilde_meet(1:end/n_s) .* int_surplus_buy(1:end/n_s));...
    sum(par.P_tr_meet(end/n_s+1:end) .* par.P_tr_tilde_meet(end/n_s+1:end) .* Pss_tilde_meet(end/n_s+1:end) .* int_surplus_buy(end/n_s+1:end))]...
    + [sum(par.P_tr_meet(1:end/n_s) .* par.P_tr_tilde_meet(1:end/n_s) .* Pss_tilde_meet(1:end/n_s) .* int_surplus_sell(1:end/n_s));...
    sum(par.P_tr_meet(end/n_s+1:end) .* par.P_tr_tilde_meet(end/n_s+1:end) .* Pss_tilde_meet(end/n_s+1:end) .* int_surplus_sell(end/n_s+1:end))]...
    - par.iota * Z * XI / w;

if par.endo_entry ==1

    fval(1:7) = [(1 / par.BETA - 1 + par.DELTA) / (1 -  par.tau_k)...
     - B * (par.P * par.q_average_prod...
     + ALPHA_0.* ALPHA * (1 - par.THETA) .* q_capital...
     + ALPHA_0.* ALPHA * par.THETA * CHI_k  * (1 - par.DELTA) .* q_money...
     + ALPHA_0.* ALPHA * par.THETA * CHI_PI .* q_profit);
     
     (par.iota / (1 - par.tau_k) - B * ALPHA * par.THETA * q_money);
     
     goods_mkt;

     ALPHA_0 - 0.5 * (1 + erf((log(q_surplus_entry) - par.c_mu) / par.c_sig / sqrt(2))); % Euler equations + goods mkt clearing    
     ];
 
else
    fval(1:7) = [(1 / par.BETA - 1 + par.DELTA) / (1 -  par.tau_k)...
     - B * (par.P * par.q_average_prod...
     + ALPHA_0.* ALPHA * (1 - par.THETA) .* q_capital...
     + ALPHA_0.* ALPHA * par.THETA * CHI_k  * (1 - par.DELTA) .* q_money...
     + ALPHA_0.* ALPHA * par.THETA * CHI_PI .* q_profit);
     
     (par.iota / (1 - par.tau_k) - B * ALPHA * par.THETA * q_money);
     
     goods_mkt;
     
     ALPHA_0 - par.ALPHA_0; % shutting down endogenous entry

    ];

end

fval = fval * 1000;



%%%% other eqm variables; only compute if eqm_all = 1
if eqm_all == 1
    H       = (B / (1 - par.ETA))^ (1 / par.ETA) * K_bar;

    int_prob_aqc    =  arrayfun(@(L, a_s, a_b, K_K_tilde) integral2(@(eps_t, eps)fun.f(eps, a_b) .* fun.f(eps_t, a_s),...
                    fun.l_bound(a_s),  min(max(L / (1 - (1 + K_K_tilde) * CHI_PI), fun.l_bound(a_s)), fun.h_bound(a_s)),...
                     @(eps_t)eps_t, @(eps_t)fun.threshold(eps_t, L, CHI_PI, CHI_PI_k, K_K_tilde),...
                    'method', int_method, 'RelTol', Rel_Tol, 'AbsTol', Abs_Tol),...
                    L_b, a_tilde_meet, a_meet, K_meet./K_tilde_meet); %note for K_K_tilde, it means buyers' capital / sellers' capital
    
    int_prob_target =  arrayfun(@(L, a_s, a_b, K_K_tilde) integral2(@(eps, eps_t)fun.f(eps, a_s) .* fun.f(eps_t, a_b),...
                    fun.l_bound(a_s), min(max(L / (1 - (1 + K_K_tilde) * CHI_PI), fun.l_bound(a_s)), fun.h_bound(a_s)),...
                    @(eps)eps, @(eps)max(fun.threshold(eps, L, CHI_PI, CHI_PI_k, K_K_tilde),eps),...
                    'method', int_method, 'RelTol', Rel_Tol, 'AbsTol', Abs_Tol),...
                    L_s, a_meet, a_tilde_meet, K_tilde_meet./K_meet); %note for K_K_tilde, it means buyers' capital / sellers' capital


    int_average_2nd_price = 2 * arrayfun(@(L, a_s, a_b) integral2(@(eps_t, eps) ((1 - par.tau_k) * B...
            * ((1 - par.THETA) * eps + par.THETA * eps_t) + 1 - par.DELTA) .* fun.f(eps_t, a_s) .* fun.f(eps, a_b),...
            fun.l_bound(a_s), fun.h_bound(a_s), @(eps_t)eps_t, fun.h_bound(a_b),'method',int_method, 'RelTol', Rel_Tol, 'AbsTol', Abs_Tol),...
            L_s, a_meet, a_tilde_meet);
    
    int_HT_credit = (1 - par.tau_k) * B * CHI_PI * (K_meet + K_tilde_meet .*...
        arrayfun(@(a_s, a_b, Z_K_tilde) integral2(@(eps_t, eps)fun.q_fun(eps_t, eps, B, CHI_q, CHI_PI, CHI_k, Z_K_tilde) .* fun.f(eps_t, a_s) .* fun.f(eps, a_b),...
        fun.l_bound(a_s), fun.h_bound(a_s), @(eps_t)eps_t, fun.h_bound(a_b), 'method', int_method, 'RelTol', Rel_Tol, 'AbsTol', Abs_Tol),...
        a_tilde_meet, a_meet, Z_meet ./ K_tilde_meet));
    
    %%%% using debt first   
    int_debt = arrayfun(@(a_s, a_b, Z_K_tilde) integral2(@(eps, eps_t) fun.debt_fun(eps, eps_t, a_s, a_b, B, CHI_q, CHI_PI, CHI_k, Z_K_tilde), fun.l_bound(a_s), fun.h_bound(a_s),...
       @(eps)eps, fun.h_bound(a_b), 'method', int_method), a_meet, a_tilde_meet, Z_tilde_meet ./ K_meet);


    %%%% coefficient of variation
    %weight     = ALPHA_0 * ALPHA;
    
    int_moment_1   = arrayfun(@(a, L, weight) (1 - weight) * exp(par.log_eps_mu + par.log_eps_sig^2 / 2) +...
        weight * integral(@(x) x .* fun.g(x, L, a), fun.l_bound(a), fun.h_bound(a)), a_meet, L_s, ALPHA_0_meet * ALPHA);

    int_moment_2   = arrayfun(@(a, L, weight) (1 - weight) * ((exp(par.log_eps_sig^2) - 1) * exp(2 * par.log_eps_mu + par.log_eps_sig^2) + (exp(par.log_eps_mu + par.log_eps_sig^2 / 2))^2)+ ...
        weight * integral(@(x) x.^2 .* fun.g(x, L, a), fun.l_bound(a), fun.h_bound(a)), a_meet, L_s, ALPHA_0_meet * ALPHA);

    int_prod_std   = sqrt(int_moment_2 - int_moment_1.^2) ./ int_moment_1; 

%     prod_std = [sum(par.P_tr_meet(1:end/n_s) .* par.P_tr_tilde_meet(1:end/n_s) .* par.Pss_meet(1:end/n_s) .* int_prod_std(1:end/n_s));...
%     sum(par.P_tr_meet(end/n_s+1:end) .* par.P_tr_tilde_meet(end/n_s+1:end) .* par.Pss_meet(end/n_s+1:end) .* int_prod_std(end/n_s+1:end))];
end


end